#################################################################################
# Calculate ACDEs for Further Analyses
#################################################################################
g.ACDE.fun <- function(data){ 
  
  ACDE <- matrix(NA, 10000, 5)
  
  for(i in 1:5){
    for(s in 1:10000){ 
      
      first <- lm(rel_deserv ~ I(liberal-i)*rel_banperm + fem + age + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=data)
      first.sim <- sim(first, 1)
      
      demediate <- coef(first.sim)[,2]%*%t(as.vector((data$liberal-i))) + coef(first.sim)[,13]%*%t(as.vector((data$liberal-i)*(data$rel_banperm))) 
      
      direct <- lm( I(rel_deserv - t(demediate)) ~ rel_banperm  + fem + age + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=data)
      direct.sim <- sim(direct, 1)
      
      ACDE[s,i] <- coef(direct.sim)[2]
      
    }
  } 
  return(ACDE)
}

###########################################################################
### Calculate ACDE for m=5 imputed data sets

ACDE.1 <- g.ACDE.fun(imp.data.1)
ACDE.2 <- g.ACDE.fun(imp.data.2)
ACDE.3 <- g.ACDE.fun(imp.data.3)
ACDE.4 <- g.ACDE.fun(imp.data.4)
ACDE.5 <- g.ACDE.fun(imp.data.5)  

c.1 <- apply(ACDE.1, 2, mean)
c.2 <- apply(ACDE.2, 2, mean)
c.3 <- apply(ACDE.3, 2, mean)
c.4 <- apply(ACDE.4, 2, mean)
c.5 <- apply(ACDE.5, 2, mean)

s.1 <- apply(ACDE.1, 2, sd)
s.2 <- apply(ACDE.2, 2, sd)
s.3 <- apply(ACDE.3, 2, sd)
s.4 <- apply(ACDE.4, 2, sd)
s.5 <- apply(ACDE.5, 2, sd)

###########################################################################
# Combine ACDE of imputed data sets

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
combined.res.des <- round(combined.res, 2)

####################################################################
g.ACDE.fun <- function(data){ 
  
  ACDE <- matrix(NA, 10000, 5)
  
  for(i in 1:5){
    for(s in 1:10000){ 
      
      first <- lm(rel_vote ~ I(liberal-i)*rel_banperm + fem + age + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=data)
      first.sim <- sim(first, 1)
      
      demediate <- coef(first.sim)[,2]%*%t(as.vector((data$liberal-i))) + coef(first.sim)[,13]%*%t(as.vector((data$liberal-i)*(data$rel_banperm))) 
      
      direct <- lm( I(rel_vote - t(demediate)) ~ rel_banperm  + fem + age + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=data)
      direct.sim <- sim(direct, 1)
      
      ACDE[s,i] <- coef(direct.sim)[2]
      
    }
  } 
  return(ACDE)
}

###########################################################################
### Calculate ACDE for m=5 imputed data sets

ACDE.1 <- g.ACDE.fun(imp.data.1)
ACDE.2 <- g.ACDE.fun(imp.data.2)
ACDE.3 <- g.ACDE.fun(imp.data.3)
ACDE.4 <- g.ACDE.fun(imp.data.4)
ACDE.5 <- g.ACDE.fun(imp.data.5)  

c.1 <- apply(ACDE.1, 2, mean)
c.2 <- apply(ACDE.2, 2, mean)
c.3 <- apply(ACDE.3, 2, mean)
c.4 <- apply(ACDE.4, 2, mean)
c.5 <- apply(ACDE.5, 2, mean)

s.1 <- apply(ACDE.1, 2, sd)
s.2 <- apply(ACDE.2, 2, sd)
s.3 <- apply(ACDE.3, 2, sd)
s.4 <- apply(ACDE.4, 2, sd)
s.5 <- apply(ACDE.5, 2, sd)

###########################################################################
# Combine ACDE of imputed data sets

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
combined.res.vote <- round(combined.res, 2)

####################################################################
g.ACDE.fun <- function(data){ 
  
  ACDE <- matrix(NA, 10000, 5)
  
  for(i in 1:5){
    for(s in 1:10000){ 
      
      first <- lm(rel_office ~ I(liberal-i)*rel_banperm + fem + age + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=data)
      first.sim <- sim(first, 1)
      
      demediate <- coef(first.sim)[,2]%*%t(as.vector((data$liberal-i))) + coef(first.sim)[,13]%*%t(as.vector((data$liberal-i)*(data$rel_banperm))) 
      
      direct <- lm( I(rel_office - t(demediate)) ~ rel_banperm  + fem + age + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=data)
      direct.sim <- sim(direct, 1)
      
      ACDE[s,i] <- coef(direct.sim)[2]
      
    }
  } 
  return(ACDE)
}

###########################################################################
### Calculate ACDE for m=5 imputed data sets

ACDE.1 <- g.ACDE.fun(imp.data.1)
ACDE.2 <- g.ACDE.fun(imp.data.2)
ACDE.3 <- g.ACDE.fun(imp.data.3)
ACDE.4 <- g.ACDE.fun(imp.data.4)
ACDE.5 <- g.ACDE.fun(imp.data.5)  

c.1 <- apply(ACDE.1, 2, mean)
c.2 <- apply(ACDE.2, 2, mean)
c.3 <- apply(ACDE.3, 2, mean)
c.4 <- apply(ACDE.4, 2, mean)
c.5 <- apply(ACDE.5, 2, mean)

s.1 <- apply(ACDE.1, 2, sd)
s.2 <- apply(ACDE.2, 2, sd)
s.3 <- apply(ACDE.3, 2, sd)
s.4 <- apply(ACDE.4, 2, sd)
s.5 <- apply(ACDE.5, 2, sd)

###########################################################################
# Combine ACDE of imputed data sets

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
combined.res.off <- round(combined.res, 2)

